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1.  Introduction 


Although  the  development  of  new  energetic  materials  (EMs)  with  decreased 
sensitivity  to  initiation  and  enhanced  explosive  power  is  of  paramount  importance 
to  the  Army,  the  discovery  of  new  EMs  that  exceed  federally  mandated  standards 
for  safety  and  stability  has  proven  difficult.  Historically,  the  search  for  new  EMs 
has  focused  on  synthesis  of  compounds  that  share  a  common  chemistry  and/or 
structural  motif  with  existing  explosives.  However,  a  novel  synthetic  approach 
based  on  the  crystal  engineering  strategy  known  as  “co-crystallization”  has  recently 
been  demonstrated  as  a  promising  route  toward  development  of  new  classes  of  EMs 
with  improved  properties.1,2  Co-crystals  are  solid  state  materials  comprising  2  (or 
more)  molecular  components  that  interact  via  noncovalent  interactions  to  form  an 
extended  network.  As  a  result  of  the  intermolecular  interactions,  the  co-crystals 
often  have  properties  that  differ  from  those  of  the  individual  components,  and  this 
feature  has  been  exploited,  with  great  success,  in  the  medicinal  community  to 
enhance  the  solubility,  stability,  and  bioavailability  of  active  phannaceutical 
ingredients.3  Given  the  positive  impact  of  co-crystallization  on  the  material 
properties  of  pharmaceuticals,  co-crystals  of  a  wide  variety  of  EMs  such  as 
2,4,6,8,10,12-hexanitrohexaaza-isowurtzitane  (CL20)  with  the  following 
compounds  have  now  been  synthesized: 

•  2,4,6-trinitrotoluene  (TNT)1 

•  1, 3 -dinitrobenzene  (DNB)4 

•  Benzotrifuroxan  (BTF)5 

•  l,3,5,7-tetranitro-l,3,5,7-tetrazacyclooctane  (HMX)6 

To  date,  much  of  the  experimental  co-crystal  research  has  focused  on  synthesis  of 
new  compounds  using  different  combinations  of  EMs.  Although  a  large  number  of 
EM  co-crystals  have  been  realized,  there  is  a  dearth  of  literature  detailing  their 
material  properties  with  most  experimental  studies  to  date  focused  on  synthesis  and 
co-crystal  structure.  There  have  also  been  experimental  analyses  of  co-crystal 
sensitivity  to  impact  and  the  relation  of  the  co-crystal’s  stability  to  that  of  the 
individual  components.6,7 

However,  other  critical  material  properties  such  as  the  response  of  EM  co-crystals 
to  increasing  temperature,  or  mechanical  defonnation,  have  not  been  determined 
experimentally,  and  only  a  limited  number  of  theoretical  studies  have  been 
reported.  Liu  et  al.  used  density  functional  theory  (DFT)  to  perfonn  a  hydrostatic 
compression  study  (0-100  GPa)  of  the  CL20-HMX  co-crystal  with  analyses  of  the 
unit  cell  volume,  band  structure,  elastic  coefficients,  and  optical  absorption 
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properties.8  Zhou  et  al.  studied  the  intermolecular  interactions  of  a  series  of  CL20 
co-crystals  using  the  Atoms  in  Molecules  approach  and  established  relations 
between  the  strength  of  intermolecular  interactions  and  the  observed  sensitivity  of 
the  co-crystals.9  Sun  et  al.10  used  the  COMPASS11  force  field  to  simulate  the 
structure,  energetic,  and  mechanical  properties  of  CL20-HMX  and  found  that  the 
predicted  bulk  modulus  of  the  co-crystal  was  smaller  than  that  of  the  pure 
components. 

Although  the  aforementioned  hydrostatic  compression  study  of  Liu  et  al.  was  based 
on  first  principles  DFT,  it  was  done  using  static  geometry  optimization  algorithms, 
which  correspond  to  zero  temperature.  It  is  well  known  that  thermal  effects, 
particularly  in  molecular  crystals,  can  often  be  pronounced.  As  an  example,  the 
experimentally  observed  a— >7  phase  transition  in  cyclotrimethylene  trinitramine 
(RDX)  does  not  occur  when  using  quantum  mechanical  (QM)  potentials  at  0  K.12 
However,  inclusion  of  thermal  effects  via  finite  temperature  quantum  molecular 
dynamics  (QMD)  simulation  does  indeed  yield  the  phase  transition  at  a  pressure  in 
excellent  agreement  with  the  experiment.13  Whereas  the  aforementioned  study  of 
Sun  et  al.10  did  rely  upon  finite  temperature  molecular  dynamics  techniques,  the 
interatomic  forces  used  to  propagate  the  atomic  degrees  of  freedom  were  obtained 
using  the  COMPASS  force  field.11  COMPASS  is  an  empirical  model,  therefore  it 
may  not  be  descriptive  of  high  temperature  and  pressure  structural  transformations 
that  may  occur  in  these  novel  systems.  Further,  as  COMPASS  is  a  nonreactive 
potential,  it  has  no  capability  to  model  thermal  or  mechanically  induced  molecular 
decompositions,  which  are  obviously  of  paramount  importance  for  explosive 
energy  release. 

In  this  work,  we  report  results  of  hydrostatic  compression  of  4  CL20-based  co¬ 
crystals  including  CL20-DNB,4  CL20-BTF,5  CL20-HMX,6  and  CL20-TNT1 
using  QMD  simulation.  We  also  include  results  for  the  epsilon  phase 
(thermodynamically  stable  at  standard  temperature  and  pressure)  of  CL20  for 
comparison.14  For  each  crystal  we  provide  the  ambient  and  high-pressure  unit  cell 
parameters,  the  bulk  modulus  (and  its  pressure  derivative),  and  shock,  particle,  and 
sound  velocity  data.  To  our  knowledge,  the  high-pressure  data  provided  in  this 
report  (with  the  exception  of  S-CL2015)  are  not  currently  available  experimentally. 

2.  Computational  Methods 

The  unit  cell  geometry  (Fig.  1)  of  S-CL2014  and  all  co-crystals 1,4-6  were  taken  from 
experiment  and  the  pressure  versus  volume  (PV)  isotherms,  at  a  temperature  of 
300  K,  for  hydrostatic  pressures  up  to  50  GPa  were  obtained  using  QMD 
simulations.  The  QMD  trajectories  (after  an  initial  geometry  optimization)  were 
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integrated  using  the  JARVIS16  software  package,  which  employs  a  leapfrog 
algorithm17  for  the  atomic  coordinates  with  external  pressure  imposed  via  a 
Berendsen18  barostat.  At  each  time  step,  the  QM  energy  and  forces  were  evaluated 
with  the  CP2K19  software  package  (executed  via  an  external  call  from  JARVIS) 
using  the  revised  Perdew-Burke-Ernzerhof20  density  functional  in  a  triple  zeta 
valence  plus  polarization  Gaussian  basis  and  the  D3  dispersion  correction  of 
Grimme  et  al.21  The  D3  correction  was  chosen  based  on  its  level  of  agreement  with 
benchmark  interaction  energies  for  gas-phase  dimers  of  several  energetic  molecules 
including  1,1 -diamino-2, 2-dinitroethylene  (FOX-7),  nitrobenzene,  and 
ethylenedinitramine,  obtained  using  coupled  cluster  theory.22  Each  simulation  was 
run  for  5  ps  (using  a  0.5-fs  time  step)  with  time  averaged  quantities  accumulated 
over  the  final  picosecond  of  the  simulation  for  all  pressures. 


E-CL20 


CL20-BTF 


CL2Q-TNT 


CL2Q-HMX 


CL2Q-DNB 


Fig.  1  Experimental  unit  cell  structures  of  E-CL20  and  co-crystals.  For  each  structure,  the 
CL20  molecules  are  red  and  the  guest  molecules  are  blue. 
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3.  Results 


3.1  PV  Data 


The  predicted  unit  cell  structures  at  zero  pressure  and  300  K  are  compared  with  the 
corresponding  experimental  values  in  Table  1.  The  agreement  between  the 
computed  and  experimental  data  for  all  systems  is  good  with  a  maximum  deviation 
of  1.8%  for  the  volume  of  the  CL20-DNB  co-crystal.  The  computed  unit  cell 
geometry  as  a  function  of  pressure  is  reported  in  Tables  2  and  3  and  the  data  indicate 
an  anisotropic  pressure  response  for  each  system.  As  an  example,  for  the 
CL20-HMX  co-crystal,  the  a  and  c  cell  vectors  show  essentially  equivalent  degrees 
of  compression  up  to  50  GPa  (-13.5%)  whereas  the  b  vector  reduces  by  almost 
20%.  Interestingly,  our  data  for  CL20-HMX  differ  from  the  theoretical  predictions 
of  Liu  et  al.8  where  it  was  reported  that  the  CL20-HMX  co-crystal  is  “much  stiffer” 
along  a  than  the  other  directions.  (As  an  aside,  it  should  be  noted  that  the  b  and  c 
axes  are  either  incorrectly  labeled  in  Fig.  2  of  Liu  et  al.8  or  the  DFT  results  obtained 
in  that  study  are  grossly  incorrect.)  The  different  observations  between  Liu  et  al.8 
and  the  present  study  could  be  attributable  to  thermal  effects,  which  were  included 
in  the  present  work  whereas  the  results  of  Liu  et  al.8  correspond  to  zero  temperature. 

Table  1  Zero-pressure  unit  cell  structure  for  S-CL20  and  co-crystals.  For  each  entry,  the  first 
number  is  the  simulated  value  (300  K),  the  second  number  is  the  experimental  value,3  and  the 
third  number  is  the  percent  difference  between  simulation  and  experiment.  Lengths  in 
angstroms,  angles  in  degrees,  volume  in  cubic  angstroms. 


System 

a 

b 

c 

a 

P 

r 

Volume 

8.807 

12.474 

13.323 

89.99 

105.36 

90.01 

1411.386 

E-CL20 

8.791 

12.481 

13.285 

90.00 

106.55 

90.00 

1397.290 

0.18 

-0.06 

0.29 

-0.01 

-1.12 

0.01 

1.01 

9.226 

11.783 

21.907 

90.00 

90.04 

89.91 

2381.539 

CL20-BTF 

9.275 

11.946 

21.577 

90.00 

90.00 

90.00 

2390.713 

-0.53 

-1.36 

1.53 

0.00 

0.04 

-0.10 

-0.38 

9.555 

13.258 

33.230 

90.00 

89.94 

90.01 

4209.337 

CL20-DNB 

9.470 

13.459 

33.620 

90.00 

90.00 

90.00 

4285.199 

0.89 

-1.49 

-1.16 

0.00 

-0.06 

0.01 

-1.77 

16.289 

9.854 

12.344 

89.99 

99.21 

90.02 

1955.885 

CL20-HMX 

16.346 

9.936 

12.142 

90.00 

99.23 

90.00 

1946.530 

-0.35 

-0.83 

1.66 

-0.01 

-0.02 

0.02 

0.48 

9.719 

19.320 

24.736 

90.00 

90.01 

90.01 

4644.368 

CL20-TNT 

9.674 

19.369 

24.690 

90.00 

90.00 

90.00 

4626.221 

0.46 

-0.25 

0.19 

0.00 

0.01 

0.02 

0.39 

a  Experimental  temperatures:  E-CL20/100K,  CL20-BTF/293K,  CL20-DNB/293K.CL20-HMX/95K, 
CL20-TNT/95K 
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Table  2  Unit  cell  geometry  (300  K)  as  a  function  of  pressure  for  E-CL20,  CL20-BTF,  and 
CL20-DNB.  CL20-BTF  simulations  were  not  convergent  at  50  GPa. 


S-CL20 


Pressure 

a 

b 

c 

a 

P 

7 

Volume 

0.5 

8.746 

12.338 

13.209 

89.98 

105.26 

90.01 

1375.141 

1 

8.697 

12.226 

13.112 

89.97 

105.11 

90.01 

1346.023 

1.5 

8.651 

12.128 

13.028 

89.95 

104.94 

90.01 

1320.622 

2 

8.620 

12.001 

12.956 

89.99 

104.68 

90.03 

1296.450 

5 

8.471 

11.621 

12.604 

89.99 

103.91 

90.03 

1204.361 

10 

8.291 

11.289 

12.255 

90.01 

103.41 

90.02 

1115.847 

20 

8.110 

10.878 

11.800 

90.02 

103.14 

89.98 

1013.662 

30 

7.963 

10.647 

11.502 

90.00 

102.90 

90.00 

950.475 

40 

7.854 

10.468 

11.273 

90.00 

102.73 

89.98 

903.978 

50 

7.758 

10.336 

11.074 

90.00 

102.54 

89.99 

866.761 

CL20- 

-BTF 

0.5 

9.125 

11.660 

21.770 

89.98 

90.06 

89.93 

2316.269 

1 

9.005 

11.568 

21.675 

90.02 

90.01 

89.96 

2257.862 

1.5 

8.932 

11.483 

21.582 

90.00 

90.05 

90.03 

2213.621 

2 

8.866 

11.408 

21.518 

90.03 

90.04 

89.93 

2176.445 

5 

8.607 

11.096 

21.111 

89.99 

90.01 

89.94 

2016.183 

10 

7.883 

11.031 

21.115 

90.25 

88.69 

88.84 

1835.101 

20 

7.542 

10.672 

20.589 

90.42 

88.65 

88.74 

1656.175 

30 

7.272 

10.474 

20.291 

90.43 

88.85 

88.79 

1544.905 

40 

6.797 

10.590 

20.174 

89.32 

90.18 

91.81 

1451.257 

50 

CL20- 

-DNB 

0.5 

9.425 

13.169 

32.958 

90.01 

89.99 

90.00 

4090.529 

1 

9.327 

13.127 

32.701 

89.99 

89.97 

90.01 

4003.731 

1.5 

9.247 

13.099 

32.433 

90.00 

90.03 

90.00 

3928.662 

2 

9.188 

13.061 

32.198 

90.02 

90.00 

89.98 

3863.821 

5 

8.929 

12.857 

31.147 

90.02 

89.99 

90.01 

3575.463 

10 

8.690 

12.595 

30.143 

90.00 

90.00 

89.99 

3299.152 

20 

8.379 

12.205 

29.030 

89.99 

90.00 

90.04 

2968.959 

30 

8.186 

11.918 

28.399 

89.96 

89.98 

89.97 

2770.411 

40 

8.036 

11.727 

27.865 

90.09 

89.94 

90.04 

2626.031 

50 

7.898 

11.563 

27.535 

89.98 

90.28 

89.95 

2514.604 
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Table  3  Unit  cell  geometry  (300  K)  as  a  function  of  pressure  for  CL20-HMX  and  CL20- 
TNT.  CL20-TNT  simulations  were  not  convergent  at  50  GPa. 


CL20-HMX 

Pressure 

a 

b 

c 

a 

P 

Y 

Volume 

0.5 

16.177 

9.763 

12.226 

89.97 

99.10 

90.00 

1906.628 

1 

16.096 

9.695 

12.133 

90.01 

99.08 

89.98 

1869.646 

1.5 

15.989 

9.631 

12.076 

90.02 

98.93 

90.00 

1837.104 

2 

15.946 

9.588 

11.966 

90.04 

99.08 

90.01 

1806.603 

5 

15.645 

9.285 

11.725 

90.05 

99.00 

89.98 

1682.302 

10 

15.305 

8.959 

11.512 

90.01 

98.86 

90.04 

1559.678 

20 

15.190 

7.891 

11.758 

89.97 

98.07 

90.12 

1395.449 

30 

14.841 

7.710 

11.510 

89.87 

97.73 

90.34 

1305.060 

40 

14.595 

7.594 

11.277 

89.75 

97.51 

90.50 

1239.080 

50 

14.079 

7.932 

10.690 

90.45 

95.40 

89.84 

1188.418 

CL20-TNT 

0.5 

9.649 

19.084 

24.567 

90.00 

89.98 

90.02 

4523.402 

1 

9.600 

18.822 

24.448 

90.00 

89.97 

90.00 

4417.430 

1.5 

9.547 

18.659 

24.326 

89.99 

89.99 

89.99 

4333.494 

2 

9.481 

18.517 

24.249 

89.98 

89.98 

89.99 

4257.111 

5 

9.274 

17.902 

23.749 

89.99 

90.01 

90.00 

3942.715 

10 

8.973 

17.265 

23.411 

90.00 

90.04 

90.00 

3626.657 

20 

8.619 

16.480 

22.988 

89.98 

89.99 

89.99 

3265.369 

30 

8.395 

16.009 

22.634 

90.00 

90.00 

90.01 

3041.775 

40 

8.242 

15.656 

22.316 

89.98 

89.99 

90.05 

2879.749 

50 

V/Vo 


Fig.  2  PV  data  for  E-CL20  and  co-crystals.  Experimental  data  taken  from  Gump 
et  al.15 

Approved  for  public  release;  distribution  is  unlimited. 

6 


Plots  of  the  computed  unit-cell  volume  as  a  function  of  pressure  are  presented  in 
Fig.  2.  As  shown,  the  pressure  response  of  the  co-crystals  is  very  similar  to  that  of 
S-CL20  and  only  show  a  slight  reduction  in  stiffness  at  high  pressure.  At  this  time, 
we  cannot  conclusively  comment  on  the  presence  of  pressure-induced  phase 
transitions  in  the  co-crystals.  The  current  data  show  no  evidence  of  high-pressure 
phase  transitions,  but  longer  time  simulations  with  larger  computational  cells  may 
show  different  behavior.  This  will  be  the  subject  of  a  future  study  using  linear 
scaling  QM  techniques23  with  first  principles  DFT  or  using  semi-empirical 
Hamiltonians  such  as  tight-binding24  to  simulate  larger  cells  for  longer  periods. 

3.2  Derived  Data 


The  bulk  modulus  and  approximations  to  the  Hugoniot  shock  and  particle  velocity 
values  were  obtained  using  the  isothennal  PV  data.  Table  4  contains  the  bulk 
modulus  (K),  and  pressure  derivative  (K’),  of  each  co-crystal. 

Table  4  Bulk  modulus  (K)  and  pressure  derivative  (K’)  of  CL20  co-crystals 


Crystal 

K(GPa) 

K’ 

CL20-BTF 

15.27 

8.45 

CL20-DNB 

15.78 

9.54 

CL20-HMX 

18.53 

8.02 

CL20-TNT 

16.98 

6.97 

E-CL20 

17.81 

7.46 

s-CL20(expt)a 

16.88 

3.82 

Experimental  data  limited  to  points  below  2  GPa  for 


comparison. 


The  moduli  were  determined  by  fitting  all  pressure  points  (up  to  2  GPa)  to  the  third 
order  Birch-Mumaghan25  equation  of  state.  With  the  exception  of  the  CL20-HMX 
co-crystal,  all  systems  show  a  reduction  in  the  value  of  the  bulk  modulus  relative 
to  the  E-CL20  reference.  The  predicted  value  of  18.53  GPa  for  the  CL20-HMX  co¬ 
crystal  is  to  be  compared  with  the  reported  value  of  8.3  GPa  reported  in  Sun  et  al. 10. 
The  value  in  Sun  et  al. 10  was  obtained  using  the  COMPASS  force  field  whereas  the 
result  reported  in  the  present  work  was  obtained  using  an  ab  initio  potential. 
Although  one  would  generally  consider  an  ab  initio  potential  to  be  more  accurate, 
our  values  may  be  suffering  from  size  effects,  which  can  lead  to  materials  that  are 
too  stiff  under  load.  Thus,  our  predicted  values  may  be  too  large.  An  experiment 
will  be  necessary  to  resolve  the  differences  between  this  work  and  the  previous 
study. 

Table  5  contains  the  predicted  sound  speed,  shock  (Us),  and  particle  (Up)  velocities 
that  were  obtained  from  the  PV  data  using  the  standard  relations.26  These  values 
represent  approximations  to  the  true  Hugoniot  loci  since  the  PV  data  in  this  study 
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correspond  to  300  K,  and  the  associated  temperature  increase  that  would  result  from 
shock  compression  was  not  incorporated.  With  that  caveat,  a  plot  of  the 
corresponding  shock  and  particle  velocities  is  given  in  Fig.  3  and,  as  shown,  all 
curves  are  essentially  linear.  However,  as  stated  previously,  using  the  current 
simulation  length  and  time  scales,  we  cannot  at  this  time  conclusively  exclude  the 
existence  of  phase  transitions  which  would  appear  as  discontinuities  in  the  Us-Up 
curves  in  Fig.  3. 

Table  5  Predicted  sound,  shock,  and  particle  velocity  values  of  co-crystals.  All  values 
obtained  using  isothermal  compression  data  from  QMD  simulation  at  300  K.  All  velocities  in 
km/s. 


CL20- 

-BTF 

CL20-DNB 

CL20-HMX 

CL20-TNT 

Pressure 

Us 

Up 

Us 

Up 

Us 

Up 

Us 

Up 

0.5 

0.08 

3.07 

0.09 

3.04 

0.08 

3.16 

0.08 

3.18 

1 

0.16 

3.16 

0.16 

3.27 

0.15 

3.38 

0.16 

3.28 

1.50 

0.23 

3.32 

0.23 

3.43 

0.21 

3.52 

0.23 

3.43 

2 

0.30 

3.47 

0.29 

3.56 

0.28 

3.62 

0.30 

3.55 

5 

0.63 

4.12 

0.63 

4.17 

0.59 

4.24 

0.63 

4.17 

10 

1.09 

4.76 

1.06 

4.92 

1.01 

4.98 

1.07 

4.90 

20 

1.78 

5.84 

1.76 

5.96 

1.70 

5.92 

1.77 

5.95 

30 

2.34 

6.66 

2.32 

6.77 

2.24 

6.73 

2.33 

6.76 

40 

2.85 

7.29 

2.80 

7.46 

2.71 

7.41 

2.83 

7.44 

50 

3.24 

8.06 

3.14 

8.00 

Sound  speed 

3.00 

3.09 

3.21 

3.10 

Fig.  3  Predicted  shock  and  particle  velocities.  Lines  have  been  added  as  guides  to  the  eye. 
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4.  Conclusion 


In  the  absence  of  experiment,  computational  modeling,  particularly  those  methods 
based  on  first  principles  quantum  mechanics,  can  provide  material  properties  that 
are  in  very  good  agreement  with  measured  data.  However,  a  true  determination  of 
the  accuracy  of  the  current  results  can  only  be  obtained  through  experimental 
validation.  Ultimately,  we  would  like  to  develop  predictive  methods  that  would 
enable  a  priori  evaluation  of  the  performance  of  a  co-crystal  based  on  the  properties 
of  the  individual  components.  The  data  provided  in  this  report  can  be  of  value  in 
the  establishment  of  these  types  of  metrics  via  comparison  of  the  co-crystal  data 
that  we  have  provided  to  that  of  the  individual  components.  However,  more  data 
for  a  much  wider  variety  of  energetic  co-crystals  will  have  to  be  computed  before 
a  general  correlation  can  be  established.  As  stated,  the  shock  and  particle  velocity 
data  in  this  work  are  approximate  since  they  were  computed  from  isothermal 
compression  values.  Precise  determination  of  true  Hugoniot  points  using  QMD 
simulation  is  currently  underway,  and  the  initial  simulations  indicate  a  significant 
increase  in  system  temperature  with  increasing  shock  pressure.  This  in  turn  drives 
chemical  reactivity  in  several  of  the  co-crystals  at  modest  pressures.  A  complete 
Hugoniot  study  will  be  the  subject  of  a  future  report. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


BTF 

benzotrifuroxan 

CL20 

2,4,6,8,10,12-hexanitrohexaaza-isowurtzitane 

DFT 

density  functional  theory 

DNB 

1 ,3-dinitrobenzene 

EM 

energetic  material 

FOX-7 

1 , 1  -diamino-2 ,2-dinitroethylene 

HMX 

1 ,3,3,7-tetranitro- 1 ,3,5,7-tetrazacyclooctane 

PV 

pressure  versus  volume 

QM 

quantum  mechanical 

QMD 

quantum  molecular  dynamics 

RDX 

cyclotrimethylene  trinitramine 

TNT 

2,4,6-trinitrotoluene 
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